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Abstract 

The construction of stochastic solutions for nonlinear partial differ- 
ential equations is a powerful method to obtain new exact results and 
to develop efficient numerical algorithms, in particular when domain 
decomposition techniques are used. 

This paper deals with the problems that arise when the nonlinear 
terms are nonpolynomial or involve derivatives. A set of equations 
of relevance for plasma physics is used as a testing ground for these 
problems. 

1 Introduction 

A stochastic solution of a linear or nonlinear partial differential equation 
is a stochastic process which, when started from a particular point in the 
domain generates after a time t a boundary measure which, integrated over 
the initial condition at t — 0, provides the solution at the point x and time 
t. For example for the heat equation 

1 d 2 

d t u(t,x) = -— u(t,x) with u(0,x) = f(x) (1) 
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the stochastic process is Brownian motion and the solution is 



u(t,x)=E x f(X t ) 



(2) 



E x meaning the expectation value, starting from x, of the process 



dX t = dB t 



(3) 



The domain here is M. x [0, t) and the expectation value in (j2J) is the inner 
product (fit, f) of the initial condition / with the measure \i t generated by 
the Brownian motion at the t— boundary. An important condition for the 
stochastic process (Brownian motion in this case) to be considered the solu- 
tion of the equation is the fact that the same process works for any initial 
condition. This should be contrasted with stochastic processes constructed 
from particular solutions. 

That the solutions of linear elliptic and parabolic equations, both with 
Cauchy and Dirichlet boundary conditions, have a probabilistic interpre- 
tation is a classical result and a standard tool in potential theory pQ [2] 
[3J. In contrast with the linear problems, explicit solutions in terms of ele- 
mentary functions or integrals for nonlinear partial differential equations are 
only known in very particular cases. Therefore the construction of solutions 
through stochastic processes, for nonlinear equations, has become an active 
field in recent years. The first stochastic solution for a nonlinear PDE was 
constructed by McKean |3] for the KPP equation. Later on, the exit mea- 
sures provided by diffusion plus branching processes [5] j6] as well as the 
stochastic representations recently constructed for the Navier-Stokes [7] [8] 
M ffD| HQ, the Vlasov-Poisson [E] [13] [IS], the Euler [TJ] and a fractional 
version of the KPP equation [16] define solution-independent processes for 
which the mean values of some functionals are solutions to these equations. 
Therefore, they are exact stochastic solutions. 

In the stochastic solutions one deals with a process that starts from the 
point where the solution is to be found, a functional being then computed 
when the process reaches the boundary. In addition to providing new exact 
results, the stochastic solutions are also, in some cases, a promising tool for 
numerical implementation. This is because stochastic simulation only grows 
with the dimension of the process, whereas a deterministic algorithm grows 
exponentially with the dimension of the space. In addition, because of the 
independence of the sample paths of the process, they are a natural choice 
for parallel and distributed computation. 
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Stochastic algorithms are also used for domain decomposition purposes 
[T7] [T%] [T9] . One decomposes the space in subdomains and then uses in 
each one a deterministic algorithm with Dirichlet boundary conditions, the 
values on the boundaries being determined by a stochastic algorithm, thus 
minimizing the time-consuming communication problem between domains. 

There are basically two methods to construct stochastic solutions. The 
first method, which will be called the McKean method, is essentially a prob- 
abilistic interpretation of the Picard series. The differential equation is writ- 
ten as an integral equation which is rearranged in a such a way that the 
coefficients of the successive terms in the Picard iteration obey a normaliza- 
tion condition. The Picard iteration is then interpreted as an evolution and 
branching process, the stochastic solution being equivalent to importance 
sampling of the normalized Picard series. The second method constructs the 
boundary measures of a measure-valued stochastic process (a superprocess) 
and obtain the solutions of the differential equation by a scaling procedure. 
For a comparison of the two methods refer to [20J . 

To extend the construction of stochastic solutions to cases more general 
than those dealt with in the past, techniques must be developed to handle 
derivatives and nonpolynomial interaction^]. Sometimes the direct handling 
of derivatives may be avoided if the derivative of the propagation kernel is 
smooth. This is the case in the configuration space Navier-Stokes equation 
[TO] , where by an integration by parts the derivative of the heat kernel is 
controlled by a majorizing kernel and absorbed in the probability measure. 
However, in general, this is not possible. 

In this paper, the construction of stochastic solutions, for differential 
equations involving derivatives and nonpolynomial interactions, will be car- 
ried out for two systems of equations which describe plasma turbulence in 
the scrape-off layer. In both cases we deal with the Cauchy problem, namely 
the equations are defined in the full space with initial conditions at t — 0. 
This is the most natural setting when the McKean approach is used. Spa- 
tial boundary conditions are easier to implement through the superprocess 
formulation, with or without a scaling limit (see [20J). 

1 Here one is concerned with the construction of stochastic solutions using McKean's 
method. Notice that the construction of Dynkin's superprocesses is also restricted to non- 
linear terms u a with < a < 2. A plausible conjecture is that, to extend the application 
of superprocesses to more general non-linear equations, one should move from processes 
on measures to processes on general distributions. 
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2 A system of scrape-off layer equations (SOLEDGE 
2D) 



The S0LEDGE-2D equations are |2T] 

d t N +-d e T + -N = Dd 2 N 
q 7] 

d t r + -(i- x )d e (^- + N)+^(r-r ) = ud 2 r (4) 

q \N J V 

where T and N are the dimensionless parallel momentum and density, (r, 9) 
are the radial and poloidal coordinates and the mask function x equals one 
in a region where an obstacle is located and zero elsewhere. 

To construct a stochastic representation for the solution one needs to iden- 
tify a stochastic process associated to the linear component (to the full linear 
component or part of it) and then, through an integral equation, construct 
the branching mechanism representing the nonlinear part. 



2.1 The linear part, \ = 

The linear part of the system for \ = is: 

d t N+-d e T = Dd 2 N 
1 

d t T + -d e N = vd 2 r T (5) 
N (0 r 6) 

Given the initial conditions at time zero I „ ' \J ) the solution of this 



system is 



r(O,r,0) 



r(^)—'0 + 1, *}(r«l (6) 



A and B being the matrices 



, , 1\ (DO 

A = [ 1 ' B = „ 
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With x = (r, 9), define a function F such that 

F(x,i) = N(x) if i = +1 

F(x,i) = T(x) if i= -1 (7) 
Associated to the equations (jSJ) there is an operator O 

(OF) (x, i) = \ {2D8 i}1 + 2v5^ x } d 2 r F (x, i) - -5 h ^d e F (x,j) (8) 
z q 

which is the generator of the stochastic process associated to the full linear 
part of the equation. 

However, for the construction of a stochastic solution to the nonlinear 
equation through a probabilistic interpretation of the integral equation, it is 
convenient to have a stochastic process that operates in a simple way on the 
arguments of the function. Therefore instead of the process associated to O, 
only the diffusion associated to the first term in (JBJ) will be used below. It 
also provides an easier handling of the dg derivative. 



2.2 The x = 1 case 

In the x — 1 case the equation (jlj) is linear 



d t N +-d T + -N = Dd 2 r N 
q r] 

d t T+-{T-T ) = vd 2 r Y (9) 
V 



the solution being 



N(t)\ _ t (_io*+«?-i) \ N(0) 

r(t) J ~ \{r(o) 

r(-IC9 e +B92_i- 



+ J dre'^-^ +a °-V ^ ) \ (10) 







with B defined before and C being the matrix 

C ( 1 
( 1 
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2.3 The nonlinear equations (% = 0): A stochastic so- 
lution 

For the nonlinear equations one writes 

1 f* 

N(t,r,6) = e tDd ?N(0,r,6)-- dTe rDd " d e T (t - T ,r,9) 

Q Jo 

-I 

Q Jo 

(11) 



I' it. r. 0) = e^r(0,r,#) - - [ dTe Tu9 "d e l^ + iv|(t-r,r,i 



Denote by £,i N ^ and £s two Brownian motions in the r— coordinate with 
diffusion coefficients \/2D and \ph). Then the equations (|TT|) may be reinter- 
preted as defining a probabilistic processes for which the expectation values 
are the functions N (t, r, 0) and Y (t, r, 9), that is 



N(t,r,6) 



E 



(t,r,0) 



1-P 



T{t,r,e) = E 



(t,r,8) 



p V / q{l-p) Jo t 

t 



drdeT (t-T,& N \6) 



(12) 



E(t,r,0) denotes the expectation value of a stochastic process started from 
(t,r,9). The processes that construct the solution at the point (t,r,8) are 
backwards-in-time processes that start from time t and propagate to time 
zero. With probability p the processes reach time zero and the contribution to 

the expectation value is ±N (o, £ t (A ° , e) (or ±r (o, & (r) , $)). With probability 

(1 — p) the process is interrupted at a time r chosen with uniform probability 
in the interval (t, 0). For the process associated to N, the process changes its 
nature, becomes a T process and picks up a factor — q ^_ p ^ ■ For the case of 
the process T, with probability |, this process either changes to a N process 
or branches into a N and a T process. In both cases it picks up a factor 

2t 

g(i-p)' 

Notice that the propagation process acts only on the r— coordinate. There- 
fore the derivative 8q, the square in T 2 and the quotient in ^ may all be 
treated as operators which are kept as labels at each branching point. When 
all the lines of the process reach time zero, the initial condition is sampled 
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at the arrival r — coordinate. This initial condition is not simply a value 
but a function of 9 (V (0,r ,6) or N (0,r Q ,6)). It implies that both the ini- 
tial condition and all its derivatives at the argument 8 must be provided. 
This initial functions are then backtracked throughout the sample lines, the 
multiplicative factors are picked up at each r interrupt and the operators 
applied whenever a labelled branching point is reached. This provides the 
contribution of each sample path to the expectation value. Figure 1 displays 
an example of a sample path, where the operators picked up along the way 
are denoted by flags. 

Notice the order of the operators at each branching point. For example, 
at the leftmost 5— labelled point the operation is 



do 



r 2 ( o, r $-\o 



Ar(o,rS 2) ,0) 



and the whole contribution of this sample path to the N— expectation value 
is 




deTlO 



times the factor (-} ttt-^i ■ 

If the initial conditions N, T and all its derivatives are bound by a 

constant M, a worst case analysis implies that almost sure convergence of 
the expectation value is guaranteed for 



t 



M < 1 



However, in practice, this condition is too severe. 



3 A two-dimensional fluid model for the scrape- 
off layer (TOKAM 2D) 

Here both the configuration and the Fourier space equations will be analyzed. 
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3.1 The configuration space equations 

A two-dimensional fluid model for the scrape-off layer based on the inter- 
change instability (TOKAM 2D [22]) is 

| n = S — {4>, n} — ane A ~^ + DA±n 

|A ± = a (1 - e A -f) + uAlcji - {0, A ± 0} - \gd x jx [i6) 

where n = -j^ is the normalized density field and 4> = y~ ^ ne norma l- 
ized electric potential. The brackets {■ • ■ } are Poisson brackets, {/, g} = 
d Xl fd X2 g — d X2 fd Xl g, with x\ = (r — a) / p s the minor radius normalized by 
the Larmor radius p\ = T e /rrii and X2 = a8/p s , a being the plasma radius. 
S is a source term. 

The first equation is rewritten as 



JL n = -ae A n + DA ± n + S - {0, n} - ae A £\ =1 ^nft (14) 

the first two (linear) terms generating the propagating part of the stochastic 
process and the others the branching. 

The second equation involves (f>(x) ,A±<f>(x) and A\cf){x). Because the 
time derivative acts on Aj_0 (x) it is natural to choose this term as the basic 
variable. Define 

* (x) = A ±( j) (x) (15) 
Then with ^ 

K{x-y) = — \og\x-y\ (16) 

(x) = (KV) (y) = j d 2 yK {x-y)V (y) (17) 
the second equation in (fT3{) becomes 



= -ae A ^ + isA ± ^ + a (l - e A ) + <re A (if + 5) * 

at 

_ ae A ^ (zlZ ( K iny _ , vl/} - gdx2 logn (18) 

where 5 stands for the kernel 5 (x — y). Notice that from the negative contri- 
bution at x = y of the term ae A K^ (in the expansion of — ae~^), a negative 
part is extracted which becomes a local dissipative term. 
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The integral version of the equations (fill) and (fTSI) is 

n(x,t) = e- t<Te V DA ^n (x, 0) 



+ jf dse- saeA e sD ^ | 5 - {if*, n} - ae A ^ ^~ p n J OM - s ) 

(19) 



tf(ar,t) = e~ taeA e tuA ^(x,0)+ [ e- saeA e svA± {a(l-e A )+ae A (K + 5)V 

Jo 

_ ae A ^ Lll! (K^y - {(if*) , *} - ^ log n 1 (x, t - s) (20) 

j=2 J - J 

The integral equations will now be given a probabilistic interpretation. In 
addition to the usual propagation and branching mechanisms, the kernel 
integrations (if*) (x) = f d 2 yK (x — y) * (y) and (d Xi K *) (x) must also be 
given a probabilistic interpretation. For this purpose two questions have to be 
dealt with. First, the kernel integrals J d 2 yK (x — y) and f d 2 yd Xi K (x — y) 
are not finite, second if (x — y) and d Xi if (x — y) are not positive definite. A 
positive function h (y) is chosen in such a way that 



N(x) = j d 2 y\K(x-y)\h(y);N 5 (x) = j d 2 y |if + 6\ h (y) ; N^x) = j d 2 y \d Xi K\ h (y) 

are Then 

\K(x-y)\h(y) |if (x - y) + 5 (x - y)\h (y) 

p{x,y) = Trr-. ;p s (x,y)- 



Pi (x, y) 



\d t K (x-y)\h(y) 



Ni(x) 

may be considered as x— dependent probability densities in y G 1R 2 . Define 

*(x,t) 



h (x) 



2 There are many functions satisfying this requirement. For example h (y) = \y\ e 
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Then, Eqs.( 1T§j) and (I2"0"j) are rewritten 

n(x,t) = e- tae \ tDA ^n{x,0)+ f dsae A e~ saeA e sDA ^ {p s , n M s , n S (x) 

Jo 

+pV n M 2>n d X2 n J d 2 y Pl (x, y) Q (y) + p%M hn d xl n J d 2 yp 2 (x, y) Q (y) 
+Px,n Pi,nM hn n (x) ^ J d 2 yp (x, y) Q (yfj J (t - s) (21) 



n (x, t) = e -' CTe e tl/A± ft (ar, 0) + J dsae A e- sae e 3 ^ 1 - |p fc M fc + p s M s J d 2 yp 5 (x, y) Q. (y) 

+Px,n Pj& M j,n ( I d2 yp ( x , y) o (y)) 

3=2 ^ ' 

Q^S^n + 5 M S°) / d2 yp* (*> y) n (v) + p 9 M ^ l °z n } (* - s ) 

(22) 

with probabilities 

Ps,n = (3 + o-e A ) 1 

Pgj. = Pg»=(3 + ^)" 1 



PE,n = ae A (3 + ae A ) 
1 



a^ 1 



Pit 



j! (e - 1) 

(2 + <7 + 2<7e A + #) _1 



= P%={2 + <j + 2<je K + g y l 

Ps = ae A (2 + a + 2ae A + g)' 1 

Pe,h = cre A (2 + a + 2ae A + #)~ 1 

p g = g (2 + a + 2ae A + g)' 1 

Hn = mhq (23) 
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and multipliers 



1 + 3/ (ae A ) 
-N^x) (1 + 3/ (ae A )) 
iV 2 (x)(l + 3/ (ae A )) 
_(_l)i( e _l) (1 + 3/ (ae A ))N(xY 

(1 - e A ) (2 + a + 2cxe A + (?) (/i (x) ae A )" 1 
-2 (2 + a + 2(xe A + (/) JVi (x) (<xe A ) _1 

2 (2 + a + 2ae A + g) N 2 (x) (ae A ) _1 
-2 (2 + a + 2(re A + JVi (x) ((xe A ) _1 d X2 log ft, (x) 
2 (2 + a + 2ae A + g) N 2 (x) (ae A ) ^ d Xl log h (x) 
iV 5 (x) (2 + (x + 2(xe A + g) (h (x) ae A ) _1 

- (2 + a + 2ae A + (/) (/i (x) ae A ) -1 

- (-I) 3 (e - 2) (/i (x) ae A )~ 1 (2 + a + 2ae A + iV (x) j (24) 

In this form the integral equations may be solved by two stochastic pro- 
cess, for n and Q, which starting from time t, propagate backwards-in-time 
as Brownian motions with diffusion coefficients D and v, respectively. These 
processes either reach time zero with probability e~ tae or branch at time 
t — s with probability density ae A e~ sae . At each branching point the appro- 
priate branching is chosen according to the probabilities p Syn , p$ n , Vbui P^,n 

for the n— process and pj., p^ n , Pbq.-> Ps, PY,,n> P g for the fi— process. With 
probability p s n the n— process samples the source S and with probability pk 
the Q— process is killed. When it is S, n or £, Q that is chosen, the branching 
has j branches according to the probability Pj >n or pj t &. For the branches of 
the n— type the argument x becomes x' = e sDA± x or x' = e suA± x, that is, it 
is propagated by the Brownian motions. For the branches of Q— type the ar- 
gument y is chosen with x— dependent probability densities p(x,y), ps (x,y) 
or pi (x,y). At each branching point the multipliers in f)24p are picked up as 
multiplicative factors. For the Q— process, the branching associated to the 
Poisson bracket term has two possibilities occurring with probability | each. 

They are either d Xi Q (x) Q (y) or Q (x) Q (y) with multipliers M^q or M^. 



M s , n 
M 1>n 



M. 



(b) 

2,n 



Mr 



M. 
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These two branchings originate from the term 



1 



d Xi (h (x) Q(x)) = Q (x) d Xi log h (x) + d Xi Q (x) 



h(x) 



This splitting is essential because the operator j^d Xi h (x) does not commute 
with the time evolution e tuA± . 

The final step in the specification of how the stochastic processes lead 
to a solution of the equations (fl9l) and ( |20l) is the handling of the deriva- 
tives d Xi n (x) and d Xi Q (x). Contrary to the case of the SOLEDGE equation 
treated in Section 2, the derivatives in this case act on the same variables as 
the stochastic process. However, d xi and d X2 commute with the time evolu- 
tion operators e* DAj - or e tuA± . Hence they may be kept as operator labels 
at the branching points and proceed with the evolution of the fields n and 
Q. When the fields finally reach time zero, the calculation is made by back- 
tracking (forward in time) through the tree the values of the fields and their 
derivatives at the final points and performing the operations at each labelled 
vertex. Notice that the operators are kept at the vertices and not carried 
along by the fields, because at non linear vertices the Leibnitz rule should 
be applied. For example, suppose that from a vertex difl (xi) at time s\ 
the field later branches into Q (#2) ^ (2/2) which then reach time zero at the 
points and x^ . Then the contribution of this derivative vertex is 



The derivatives are computed at the arrival points, whereas the terms h (x) 
and log h (x) in the multipliers (which do not commute with time evolution) 
are computed at the branching points. That this is the correct procedure is 
easily understood by recursively iterating the equations. Figure 2 shows how 
a sampled path of these processes looks like. 

The contribution of this sample path to the n— expectation value is 



M i>n (x) M jin (x) {djUQ (23) diQ (z 4 ) + didjn (z 3 ) (24)} Mfih (yi) Q (zi) d k Q (z 2 ) 



Some remarks on this construction: 

(i) The normalization of the kernels by the h (x) function converts the 
computation of the integral into a probabilistic sampling. At the same time 
the kernel itself absorbs one of the derivatives, that otherwise would act on 
the fields. 
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(ii) The exponential was expanded into a power series leading to 
multiple branchings with probabilities pj^ n and Pj^ci- By contrast the other 
nonpolynomial term logn was kept as an operator. The reason is that the 
coefficients in the series expansion of the log do not lead to a probability 
distribution. 

(iii) The complexity of the resulting processes, reflects the complexity of 
the original equations f[T3")) . Nevertheless what this construction shows is that 
clearly defined techniques may be used to construct stochastic solutions even 
when derivatives and complex non- linear terms are present. 

To insure convergence of the calculation of the solution by the stochastic 
processes, a bound M must be put on the multipliers (which depend on the 
normalization factors, hence on the h(x) function), on the initial conditions 

a(l-e A ) 

n (x) and fl (x) and their derivatives, on the source term S and on K h , < ' . 
From the branching process one estimates the probability of a tree of n 
branches, therefore 

1 



M < 



1 - e- 



e A t 



3.2 The Fourier-transformed TOKAM 2D 

In the Fourier transformed equation, the derivatives become simpler multi- 
plicative factors. However the nonlinear terms become more complex. Let 

/ = log n 

Then, the first equation in (fTBl) becomes 

^f = S'- {0, /} - ae K -t> + DA ± f + D \Vf\ 2 (25) 
where the source term has been redefined as S' — -. Let the Fourier trans- 

n 

forms of / and (f) be 

F(k,t) = -L f d 2 xf(x,t)e ik -* 



$(k,t) = — I d 2 x<j)(x,t)e ik - x 

27T J 
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Then, the Fourier transform of the equations (1131) is 

(k, t) = -D \k\ 2 F (k, t) + S' (k, t) + ^J d 2 i (h - 6) 6 (k - F (0 
-F(k 

-ae A V 1 ' f $* w (A;) — — / d 2 £ {(Jfci - £i) Ci 

+ (^-6)6)^-0^(0 
JU(M) = -HA:| 2 $(A;,t) + ^y'0(fc 1 -6)6<f(^-O$(O{IO 2 -l^-O 2 } 

(*) - ^ (*) - g -g2>- (*)} (26) 

where S' is the Fourier transform of S' and $* n denotes the convolution 
power. 

/n / n \ n 

i=i V «=i / p=i 

To control the growth of the functional obtained by the stochastic process 
one divides F (k, t) and $ (k, t) by a majorizing kernel 7 (A;) 

x(M) = « 

7 W 

C(M) - ^ P7) 
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and writes integral equations for X (k, t) 

X (k,t) = e-^l fc l 2 ^ (A:, 0) -h^'rf^ |/c| 2 e ^^l fc l 2 - ^ rf 2 ^ ( A;, ^) J> 
^{(k-U-s)x(U-s)-^x(k-U-s)C(U-s) 

+Pb^2Pn / d 2 £l ■ ■ ■ d 2 £n-lM n (k, £1, • • • , Cn-l) V (k, Cl, • • • , £n-l) 

n 



n-1 



X 



i=l 



+Pc / (*, 



+p s M s (fc) 



S'(k,t-s) 
l{k) 



X(k-£,t-s)x(£,t-s) 

(28) 



with probabilities 



Pa 
Pn 



(l + D + 27r + 27rae A ) _1 
27r ( je A+1 (1 + D + 2tt + 27r<7e A ) _1 
D(l + D + 27T + 27r ( Te A )~ 1 
2tt (1 + D + 2n + 2nae A y 1 
(eniy 1 

7(fc-Er = " 1 1 6)7(ei)---7(en-i) 

7 *n 



(29) 
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and multipliers 

M n (fc) 
M a (fc) 



(fei - 6) 67 * 7 (fc) 

7 (£;)£> |fc| 2 7r 

(fci - 6) 67 * 



(1 + £> + 2tt + 27T(7e A ) 



(k) 



(1 + D + 2tt + 2irae K ) 

(fc2 - 6)67 ? (A:) (1 + D + 2vr + 2^) 
(-^e^(k) Al + D + 2n + 2naeA) 



(2ix) 1l 1 (k)D\k\ 
1 '\ + D + 2tt + 27nxe A ) 



2ttD|A;| 



(30) 



Notice that the term n = in the sum is pbPoM 5 (k 2 ) which, when chosen, 
kills the contribution of the corresponding sample path for k 7^ 0. 
For q (k, t) 

C(k,t) = e—l fc l 2 *C(A:,0) +^rf^|A;| 2 e—l fc l 2s y rf 2 ^ (A:, P 

x((k-Z,t-s)((Z,t-s)+p' b Y,Pnj d%---d%_ 1 M: i (k^ 1 , 

X P (k^ lr -- ,C n -l)C [k-^2Zi,t-s \ ■■■C(Ca-l,t- S) 



i=l 



+ ft M x (U)x(M-s)} 
with probabilities and multipliers 







Pa = 


\2n 




(± 


P'b = 


\2tt 




(J- 


Px = 


\2ir 


Pn = 


(en!) 



1 1_ 

2tt 



7(fc-Er = 76)7te)---7(en-i) 

7 *n 



(32) 
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M , = (* 1 -{ 1 )6{ltl'-|*-tl'}7.7W fj_ 

7 (fc)z/|fc| 4 \2n 



;-l) n e A 7 w (fc) /J_ 
(2 7 r) n - 1 7 (fc)^|fc| 4 V 2 ^ 



2vr(e A - 1) / 1 
7 (/c)i/|fc| \ 27r 

^x(*) = 4^(^ + 2 + ^ ( 33 ) 



The probabilistic interpretation of these equations is similar to the previ- 
ous cases. There are two stochastic processes that started at t, propagate 
backwards in time as Brownian motions with coefficients D\k\ 2 and v |fc| 2 . 
With probabilities e _D ' fc ' * or e~ v \ k \ 1 the processes either reach t — without 
branching or branch with probability densities D |/c| 2 e _D ' fc ' s or v |/c| 2 e~ v \ k \ s 
at time t — s. The branch that they follow afterwards is controlled by the 
probabilities in (129]) and (!32l) . Because products become convolutions under 
the Fourier transform, probability densities are needed to decide the new 
momenta of the branches. The main role of the majorizing kernel 7 (k) is to 
normalize these probability densities in momentum space. At each branching 
point the process picks up the multipliers in fl30|) and fl33|) . the contribution 
of each sample path to the the expectation values of x (k, t) and ( (k, t) be- 
ing the initial values of the corresponding fields x (k, 0) or £ (k, 0) when the 
processes reach time zero multiplied by the multipliers picked up along the 
path. In this case because there are no operator labels, there is no need to 
backtrack forward in time as in the configuration space equation. The price 
one pays is a more complex branching structure at each vertex. 

A strict condition for the convergence of the process is obtained by a 
bound M < 1 on all the factors that intervene in the final calculation of each 
path contribution. 



4 Remarks and conclusions 

1) The scrape-off layer equations provide a good testing ground for the con- 
struction of stochastic solutions when the equations contain nonpolynomial 
terms and derivatives. Therefore, the techniques developed here considerably 
extend the range of equations to which the stochastic solution technique may 
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be applied. 

2) To obtain good accuracy in the stochastic solutions, many sample paths 
should be computed for each starting configuration. How many samples are 
needed may be estimated by large deviation techniques (23]. The sample in- 
dependence character of the stochastic calculation reduces the severity of this 
problem when using parallel computing. Nevertheless it should be pointed 
out that, when solutions are desired over a large domain, the stochastic 
method is computationally competitive only when applied together with the 
domain decomposition method [17] [18] [19]. Then a considerable improve- 
ment is obtained. However, if local solutions on configuration or Fourier 
space are desired, the stochastic method is quite appropriate. For example 
a study of the time evolution of a few high Fourier modes gives informa- 
tion on the turbulence spectrum that only a very fine grid and an expensive 
calculation would provide with a global deterministic algorithm. 
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Figure 1: A sample path of the N — F stochastic process 
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Figure 2: A sample path of the n— process 
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